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In  studies  of  low-frequency  reverberation  within  the  marine  environment,  a  central 
concern  is  the  relationship  between  reverberation  events  and  morphological  fea¬ 
tures  of  the  seafloor.  A  time-domain  migration  algorithm  for  the  reverberation  in¬ 
tensity  field  is  developed  that  produces  scattering  coefficient  maps  coregistered  with 
a  bathymetry  database.  The  algorithm  is  tailored  to  broadband  transient  sources 
with  good  range  resolution,  and  was  developed  to  analyze  an  extensive  set  of  rever¬ 
beration  records  from  a  200-255  Hz  source  collected  on  the  flanks  of  the  Mid- At¬ 
lantic  ridge.  The  precise,  sample-by-sample,  tracking  of  wavefronts  across  elements 
of  the  bathymetry  database  that  forms  the  foundation  of  the  algorithms  implementa¬ 
tion,  results  in  reverberation  maps  that  show  a  clear  and  detailed  correlation  be¬ 
tween  scattering  and  morphology,  with  narrow  scarp  slopes  consistently  highlighted. 
Environmentally  induced  asymmetries  in  transmission  loss  and  incidence  angle  are 
exploited  to  break  the  inherent  left-right  ambiguity  of  the  receiver  array.  Iterative 
migration,  assuming  a  dominant  dependence  of  backscatter  on  grazing  angle,  pro¬ 
duces  images,  even  from  individual  records,  that  show  good  ambiguity  resolution. 
Results  from  multiple  records  corroborate  the  effectiveness  of  the  ambiguity  resolu¬ 
tion  and  demonstrate  the  stability  of  the  scattering  coefficient  estimates  and  the 
acoustic  system. 

PACS  numbers:  43.30.Gv  43.30.  Vh  43.30.Pc 

INTRODUCTION 

In  July  of  1993,  a  coordinated  set  of  low  fre¬ 
quency  acoustics  experiments  were  conducted  on 
the  western  flank  of  the  Mid-Atlantic  Ridge  within 
a  corridor  that  lay  just  north  of  the  Kane  fracture 
zone.  These  experiments  were  designed  to  investi¬ 
gate  and  quantify  the  relationship  between  acoustic 
reverberation  and  seafloor  morphology  and  consti¬ 
tuted  the  Main  Acoustics  Experiment  (MAE)  of  a 
Special  Research  Program  (SRP)1  on  bottom  rever¬ 
beration  sponsored  by  the  Office  of  Naval  Research. 

Included  within  the  experimental  suite  were 
monostatic  and  bistatic  reverberation  measurements 
by  the  research  vessels  Cory  Chouest  and  Alliance, 
and  near  bottom  measurements  of  low  grazing  angle 
interactions  using  vertical  arrays  of  hydrophones. 

One  element  of  the  reverberation  experiments  that 
makes  them  unique  is  the  availability  of  detailed 
bathymetric  and  geophysical  data  for  the  area,  col¬ 
lected  on  a  pair  of  SRP  sponsored  large-scale  and 


small  scale  geophysics  cruises.  The  bathymetric  da¬ 
tabase  includes  a  map  of  the  entire  SRP  experimen¬ 
tal  corridor  gridded  at  200  m  and  finer  scale  sur¬ 
veys  of  selected  areas  with  resolution  on  the  order 
of  10  m  or  better .  It  is  the  availability  of  these  data¬ 
bases  that  permits  the  investigation,  in  detail,  of  the 
correspondence  between  reverberation  returns  and 
bathymetry. 

The  process  of  mapping  reverberation  returns 
back  onto  the  seafloor  scattering  sites  that  produced 
them  can  be  referred  to  as  charting2'4.  Here  we  pre¬ 
fer  to  use  the  term  migration,  because  of  the  strong 
similarity  of  the  methods  employed  here  to  the 
Kirchoff  migration  method  used  in  seismic  reflec¬ 
tion  processing.  The  principal  difference  between 
the  two  is  that  the  reverberation  data  -  beamformed, 
time-series  of  acoustic  intensity  -  are  assumed  to 
be  the  result  of  incoherent  scattering  rather  than  co¬ 
herent  reflections. 

In  outline,  migration  is  accomplished  in  two 
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stages,  first  the  seafloor  is  mapped  into  the  time- 
beam  coordinates  of  the  data,  then  the  reverbera¬ 
tion  signal  from  a  given  beam  and  small  time  inter¬ 
val  is  distributed  over  the  corresponding  ensonified 
area  of  the  seafloor.  Migration  is  complicated  by 
the  fact  that  data  were  recorded  by  a  horizontal  line 
array  (HLA)  and  each  beam  is  associated  with  a  pair 
of  directions  that  are  orientated  symmetrically  with 
respect  to  the  array  axis,  thus  a  reverberation  return 
could  originate  from  one  of  a  pair  of  distinct 
ensonified  areas.  Simple  mapping  schemes  that  ig¬ 
nore  variations  in  seafloor  depth  produce  reverbera¬ 
tion  maps  that  are  perfectly  symmetric  with  respect 
to  the  array  axis  and  provide  no  means  of  resolving, 
from  a  single  observation,  the  inherent  "left-right 
ambiguity"  introduced  by  a  HLA.  Accounting  for 
the  bathymetry  introduces  "environmental  symme¬ 
try  breaking"3,4  -  bathymetry  induced  variations  in 
transmission  loss,  acoustic  shadowing  -  which  can 
be  exploited  by  algorithms,  such  as  the  one  presented 
here,  to  reduce  substantially  the  left-right  ambigu¬ 
ity,  even  for  a  single  ping. 

An  additional  problem  with  simpler  schemes  is 
•  that  the  migrated  images  are  in  the  time-beam  coor¬ 
dinates  of  individual  transmissions  rather  than  in  glo¬ 
bal  coordinates  ef-the  bathymetry,  complicating  the 
comparison  of  multiple  pings  and  studies  of  the  re¬ 
lationship  between  scattering  and  bathymetry.  The 
current  migration  algorithm  avoids  this  problem  by 
being  formulated  directly  in  terms  of  the  backscat- 
ter  strength  of  individual  elements  of  the  bathym¬ 
etry  grid. 

The  algorithm  is  based  on  a  time-domain  for¬ 
mulation  of  the  incoherent  scattering  process  which 
relates  the  scattering  strength  to  the  expected  acous¬ 
tic  intensity  via  a  large  sparse  matrix.  This  formula¬ 
tion  means  that  we  could  treat  the  migration  as  a 
linear  inverse  problem  and  invert  multiple  pings  si¬ 
multaneously  for  scattering  strength.  Combining 
pings  with  diverse  look  angles  would  eliminate  left- 
right  ambiguity  from  the  solutions.  A  linear  inverse 
approach  would  also  allow  regularization  constraints 
to  be  employed  on  the  solutions;  an  approach  we 
pursue  elsewhere5. 

In  this  paper,  we  solve  the  scattering  equations 
approximately  using  a  version  of  iterative 
backprojection  under  the  assumption  that  backscat- 
ter  strength  is  spatially  isotropic  and  a  function  only 


of  the  local  ensonification  angle.  The  advantage  of 
using  iterative  backprojection  is  that  it  is 
computationally  fast  and  guaranteed  to  produce  mi¬ 
grated  images  that  satisfy  the  data  exactly  while  only 
weakly  enforcing  assumptions  about  the  behavior 
of  scattering  strength.  Iterative  backprojection  could 
potentially  be  used  as  a  basis  for  producing  well 
resolved,  migrated  images  in  near  real  time.  We  use 
it  here  to  rapidly  check  the  consistency  of  multiple 
images  from  successive  reverberation  records  of  the 
MAE  and  to  produce  preliminary  estimates  of  back- 
scattering  strength  as  a  function  of  grazing  angle. 
The  linear  inverse  approach  outlined  above  will  only 
be  effective  if  the  individual  reverberation  images 
are  basically  consistent.  Uncorrected  errors  in,  for 
example,  array  heading  or  source  location  would 
substantially  degrade  the  result  of  inverting  multiple 
pings. 

We  concentrate  on  the  analysis  of  a  subset  of 
the  data,  namely  half  convergence  zone  (1/2  CZ), 
monostatic  reverberations  from  broadband,  200-255 
Hz,  linearly  frequency  modulated  (LFM)  transmis¬ 
sions  recorded  by  the  RA/  Cory  Chouest .  We  chose 
this  subset  because  it  has  the  best  range  resolution 
and  simplest  propagation  paths,  and  thus  affords  the 
best  opportunity  of  investigating,  quantitatively,  the 
relationship  between  backscattering  strength  and 
bathymetry.  Previous  examinations  of  the  SRP 
acoustics  data  have,  for  the  most  part,  concentrated 
on  the  longer  range,  continuous  wave  (cw)  data,  and 
have  established  that  there  is  abroad  correspondence 
between  high  amplitude  reverberations  and  well 
ensonified  features  at  1/2  and  1 1/2  CZ,  with  the  high¬ 
est  returns  from  backfacing  ridges.  However,  the 
detailed  investigation  and  quantification  of  this  re¬ 
lationship  is  as  yet  incomplete;  this  paper  represents 
one  step  in  that  process. 

I.  OVERVIEW  OF  DATA  AND  EQUIPMENT 

The  MAE  was  conducted  within  an  approxi¬ 
mately  4°  by  2°  corridor  on  the  flanks  of  the  Mid- 
Atlantic  ridge  extending  from  45°  W  to  49°  W  and 
from  25°  30’  N  to  27°  30'  N.  The  corridor  lies  within 
a  larger  area  designated  as  an  ONR  natural  labora¬ 
tory6,  which  straddles  the  Mid-Atlantic  ridge  itself 
and  is  bounded  to  the  south  by  the  Kane  fracture 
zone.  The  seafloor  fabric  within  the  corridor  is  domi¬ 
nated  by  lightly  sedimented,  ridge  parallel  abyssal 
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FIG.  1.  Bathymetry  and  a  portion  of  FW  Cory  Chouesl  track  lines(heavy  black  lines)  in  the  vicinity  of  target  B',  an  elongated 
ridge  at  the  western  edge  of  the  Fig.  In  this  paper  we  analyze  reverberation  data  from  Pings  at  the  start  of  run  5a  and  6  of  the 
MAE  (heavy  white  circles). 

hills,  but  there  are  also  large  sedimented  ponds  and  target.  We  have  chosen  the  first  intersection  here 
the  area  is  dissected  by  deep  corridors  that  are  the  because  the  associated  transmissions  provide  better 
fossilized,  off-axis  expression  of  spreading  segment  ensonification  of  the  front  slope  of  B'  at  1/2  CZ. 
boundaries  at  the  ridge  axis.  The  feature  whose  back-  The  basic  features  of  the  RN  Cory  Chouest 
scatter  we  will  examine  in  detail  here  is  a  -30  km  source  and  receiver  arrays  are  shown  in  Fig.  2.  The 
long  abyssal  hill,  designated  B\  that  lies  at  the  west-  vertical  source  array  consisted  of  10  elements  spaced 
em  end  of  one  of  the  segment  boundary  corridors.  at  2.29  m  and  centered  at  101  m,  while  the  horizon- 

In  particular  we  will  concentrate  on  a  set  of  1/2  CZ,  tal  receiving  array  consisted  of  128  hydrophone 

monostatic  reverberation  returns  recorded  by  the  RJ  groups  spaced  at  2.5  m.  Beamforming  of  the  rever¬ 
ie  Cory  Chouest  that  span  the  intersection  of  a  pair  beralion  data  was  performed  aboard  ship  using  a 

of  the  experimental  runs,  run  5a  and  run  6,  Fig.  1  delay  and  sum  beamformer  with  a  Hamming  win- 

and  Table  I.  This  intersection  lies  slightly  northwest  dow,  resulting  in  a  nominal  broadside  beamwidth 
of  a  large  set  of  intersecting  tracks  that  constitutes  of  ~  1 .3°,  measured  at  the  -3  dB  point.  The  presence 

the  focus  of  an  extensive  set  of  bistatic  and  of  dead  or  improperly  gained  phones  in  the  HLA 

monostatic  experiments  that  had  B’  as  the  primary  resulted  in  nearly  uniform  sidelobe  level  of  approxi- 
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Run  Segment  ..  Ping  .  Location  (Deg)  :  Location  (UTM  km)  .  HLA  Array 


Lat 

Lon 

Lat 

Lon. 

Heading 

5a 

415 

1313 

26°  38.6*  N 

-47°  48.6’ W 

2950.1 

220.3 

168.6 

5a 

417 

1318 

26°  37.4'  N 

-47°  48.1’ W 

2947.8 

221.1 

166.1 

5a 

420 

1327 

26°  35.7'  N 

-47°  47.5'  W 

2944.7 

221.9 

169.5 

5a 

421 

1331 

26°  35.2’  N 

-47°  47.3' W 

2943.7 

222.3 

169.1 

5a 

423 

1336 

26°  34.0'  N 

-47°  46.9’ W 

2941.4 

222.9 

167.6 

6 

547 

1683 

26°  38.7’  N 

-47°  47.3' W 

2950.2 

222.5 

197.3 

6 

548 

1687 

26°  38.1’ N 

-47°  47.4'  W 

2949.1 

222.2 

195.9 

6 

550 

1692 

26°  36.9’  N 

-47°  47.6' W 

2946.9 

221.8 

194.8 

6 

551 

1696 

26°  36.3'  N 

-47°  47.7’ W 

2945.7 

221.6 

194.6 

6 

554 

1705 

26°  34.5'  N 

-47°  48.1' W 

2942.4 

220.9 

195.3 

Table  I:  List  of  reverberation  records,  pings,  analyzed  in  this  paper 


mately  -30  dB'. 

The  schedule  of  source  transmissions  for  the 
experiment  was  conceptually  arranged  in  a  hierar¬ 
chical  fashion  with  the  bottom  two  levels  in  the  hi¬ 
erarchy  being  transmission  segments  and  individual 
source  wavetrains  or  pings.  Segments  were  num¬ 
bered  consecutively  and  each  one  lasted  12  minutes. 
Ping  numbers  were  also  numbered  consecutively 
with  typically  6  source  wavetrains  being  transmit¬ 
ted  during  each  segment.  For  the  portions  of  runs 
5a  and  6  examined  here,  2  out  of  every  3  segments 
were  used  for  R/V  Cory  Chouest  transmissions  and 
the  remaining  segment  for  Alliance  transmissions. 
One  of  the  standard,  usually  the  first,  transmissions 
within  a  Cory  Chouest  segment  was  a  5  s.  duration 
linearly  frequency  modulated  (LFM)  signal  chirped 
over  the  band  200-255  Hz.  It  had  the  largest  band¬ 
width  of  any  of  the  source  wavetrains  and  hence  the 
best  range  resolution  after  matched  filtering,  Ar  = 
13.6  m;  it  is  thus  the  best  choice  for  examining,  in 
detail,  the  spatial  structure  of  the  scattering  process. 

The  structure  of  the  LFM  source  wavefield  out 
to  beyond  a  1/2  CZ  is  shown  in  Fig.  3.  The  broad¬ 
band  transmission  loss  was  calculated  by  integrat¬ 
ing  across  the  band,  the  transmission  losses  predicted 
at  individual  frequencies  by  a  wide-angle  parabolic 
equation  (PE)  code7.  At  times  during  the  experiment 
the  source  array  was  steered  down  at  various  angles 
with  respect  to  the  surface,  but  for  the  transmissions 
of  interest  and  the  calculation,  the  steering  angle 
was  set  at  0°.  The  sound  speed  profile,  displayed  in 
Fig.  3a,  was  essentially  the  same  throughout  the 
experiment.  The  center  of  the  source  array  at  1 8 1  m 
was  in  the  upper  part  of  the  waveguide  with  the  con- 


FIG.  2.  Sketch  of  the  source  and  receiver  array  of  the  R/V 
Cory  Chouest.  The  10  element  source  array  had  a  spacing  of 
2.29  m  while  the  126  element  receiver  array  had  a  spacing  of 
2.5  m. 

jugate  depth  being  at  3800  m  and  the  turning  range, 
1/2  CZ  distance,  being  around  33  km.  The  core  of 
the  main  acoustic  beam  contains  a  pair  of  high  am¬ 
plitude  fringes  which  result  from  the  interference 
between  energy  that  was  up  and  downgoing  at  the 
source  level.  Behind  the  core,  the  two  components 
are  sufficiently  separated  in  time,  up  to  ~30  ms,  that 
there  is  no  significant  interference  effect  on  the 
amplitude.  Furthermore  in  this  region  the  ray  propa¬ 
gation  angles  diverge  by  about  2°  (Fig.  4).  Beyond 
the  turning  range,  most  of  the  structure  in  the  main 
beam  is  due  to  caustics  in  either  the  up  or  downgoing 
source  field.  The  source  field  also  has  a  number  of 
sidelobes  that  intersect  the  seafloor  at  ranges  less 
than  15  km. 
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Range  (km) 


FIG.  3.  (a)  Representative  up  (blue)  and  downgoing  (red)  rays  paths  within  the  main  acoustics  beam  from  a  point  source  at  181 
m,  the  center  depth  of  the  vertical  source  array.  The  reference  sound  speed  (black)  has  a  surface  velocity  of  1540  m/s  and  a 
minimum  velocity  of  1494  m/s  at  1 . 1  km.  The  downgoing  ray  paths  are  simple  but  caustics  form  beyond  the  1/2  CZ  at  30-35  km. 
(b)  PE  calculation  of  the  tranmission  loss  for  a  broadband,  200-255  Hz,  source  from  the  10  element  source  array.  Transmission 
loss  variations  near  the  axis  of  the  downgoing  main  beam  are  the  result  of  interference  of  the  downgoing  and  surface  refracted 
energy.  Transmission  loss  is  more  variable  beyond  a  1/2  CZ  due  to  the  presence  of  caustics. 
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II.  MIGRATION  ALGORITHM 

The  goal  of  the  migration  algorithm  is  to  pro¬ 
duce  a  map  of  scattering  strength  registered  on  the 
•  same  grid  as  the  bathymetry  data.  The  basic  obser¬ 
vational  data  are  the  intensity  time  series  I(<pb,t)  as¬ 
sociated  with  beam  angle  (pbthat  result  from  squar¬ 
ing  the  matched  filtered  and  beamformed  pressure 
series.  The  relationship  between  the  expected  inten¬ 
sity  of  the  backscattered  field  and  the  scattering  co¬ 
efficient  is  taken  as  the  following  integral  over  the 
seafloor 

7(9  *t)  =JfdA  y  m(x,Q,§)  g{x,  cp  b,<?J-TrTs) 

(1) 

where 

y = ym/yXis^splT  (2) 

Tr,s  are  the  ray  theoretical  travel  times  to  and  from 
the  scattering  location,  yis  the  product  of  the  trans¬ 
mission  losses,  vFs>r,  to  and  from  the  seafloor;  the 
gains,  ys,  yr,  and  ymf,  associated  with  the  source  ar- 
;  ray,  receiver  array  and  matched  filtering;  and  the 
source  level,  a  product  of  the  transmission  level,  pr_ 
and  signal  duration,  T.  Since  we  are  dealing  with 
broadband  transient  signals,  gains  for  individual 
'  components  of  the  system  are  defined  in  terms  of 
the  ratio  of  integrals  over  the  input  and  output  pres¬ 
sure  fields8  - 

y=[  dtp20{t)/(  'dtpftt) 

Jo  Jo  (3) 

The  scattering  coefficient,  m(x,6,<J>)  in  Eq.  (1), 
is  a  function  of  position,  x,  and  the  monostatic  inci- 
dence/backscatter  direction  at  the  seafloor  (0,<J)).  The 
quantity  g(x,  cpb,  <pa,  t-Tr  -Ts)  can  be  termed  the 
scattering  function  for  the  system,  it  depends  on  the 
details  of  the  scattering  process,  as  well  as  quadrati- 
cally  on  the  effective  source  function  at  the  seafloor 
and  quadratically  on  the  time  domain  response  of 
the  beamformer  for  a  given  steering  angle  cpt,  and 
arrival  angle  cpa  at  the  receiver  array.  The  partition 
of  scattering  between  the  scattering  coefficient  and 
the  scattering  function  is  defined  by  normalizing  g(t) 
so  that  its  time  integral  is  unity. 

Eq.  (1)  is  a  time  domain  variant  of  the  standard 
equation  relating  the  expected  reverberation  inten¬ 
sity  to  scattering  coefficients  for  cw  signals8*9.  The 
additional  element  in  the  time  domain  formulation 


Range  (km) 


FIG.  4.  (a)  Grazing  angle  of  rays  at  a  depth  of  3.5  km.  Arriv¬ 
als  are  for  rays  from  the  top,  middle  and  bottom  elements  of 
the  source  array.  There  is  approximately  a  2°  difference  in  ray 
angle  between  surface  reflected/refracted  rays,  upper  cluster, 
and  rays  leaving  the  source  travelling  downwards.  At  this 
depth,  the  penumbra  of  the  main  beam  starts  at  25  km  and 
extends  to  about  28  km.  (b)  Difference  in  travel  times  at  3.5 
km  for  up  and  downgoing  rays  at  the  source. 

is  the  introduction  of  g,  which  can  be  regarded  as 
representing  the  resolution  function  of  the  acoustic 
imaging  system.  Analytic  expressions  or  empirical 
forms  could  be  found  for  g  (and  also  m),  given  ex¬ 
plicit  assumptions  about  the  nature  of  the  surface, 
for  example  small  amplitude  roughness8.  In  prac¬ 
tice,  we  assume  that  the  details  of  the  time  response 
will  not  be  important  provided  that  the  resulting 
scattering  coefficients  estimates  represent  averages 
over  a  sufficiently  large  area.  We  thus  use  a  simple 
factored  boxcar  form  for  g 

£(9i/9«/0=£*(0  &(9*<p.) 

=  B(f/UB((< p„-cpfl)/Acp)  (4) 
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where  B(t).is  a  unit  boxcar  centered  at  zero  and  width 
one,  and  A9  is  the  width  of  the  beam.  In  the  limit 
tw  =  0,  the  time  dependence,  g*(t),  reduces  to  a  8- 
function.  A  potentially  more  accurate  representation 
of  the  angular  dependence  would  be  to  use  the  re¬ 
sponse  function  of  the  beamformer  at  a  representa¬ 
tive  frequency  rather  than  a  boxcar . 

We  can  derive  a  sonar  equation  from  Eq.  (1)  by 
substituting  for  the  scattering  function  from  Eq.  (4). 
With  g*(t)  set  equal  to  a  8-function,  integrating  from 
t,  to  t,+T4  yields 

p\TB  =  %% . ymfysyr  p\T.  m.A  (5) 

where 

1  r1  i+T« 

P^tJ  J(<M  (6) 

and  A  is  the  ensonified  area  corresponding  to  the 
integration  limits  in  time  and  beam  angle,  and  it  is 
assumed  that  transmission  losses  etc.  are  constant 
over  A.  This  is  equivalent  to  the  following  sonar 
equation  for  the  reverberation  level  R  =  201og  pr 

R  =  -TLr  -  TLS  +  G  +  M 

+  101og(A)  +  S  -  lOlog  (T,  /  T)  (7) 

where  TLr,s  are  the  transmission  losses,  G  is  the 
combined  system  gains,  M  the  scattering  strength 
and  S  the  source  strength,  all  defined  in  the  obvious 
way  from  Eq.  (5).  The  final  term  corrects  for  differ¬ 
ences  between  the  averaging  time  and  the  signal 
duration.  Frequently  for  cw  signals  the  averaging 
time  is  chosen  equal  to  the  pulse  duration  and  there 
is  no  correction.  However,  for  the  LFM  pulse,  choos¬ 
ing  an  averaging  time  equal  to  the  width  of  the 
matched  filtered  pulse,  18  ms,  corresponds  to  a  cor¬ 
rection  of  24  dB.  A  problem  with  choosing  a  short 
averaging  time  is  that  it  is  then  hard  to  justify  the 
use  of  a  8-function  for  the  scattering  function,  since 
effectively  the  convolutional  nature  of  the  reverbera¬ 
tion  process  is  being  ignored. 

The  above  result  is  most  naturally  viewed  as 
producing  a  scattering  strength  map  in  the  local 
time-angle  coordinates  of  the  beamformer:  the 
scattering  strength  estimate  obtained  from  Eq.  (7) 
is  assigned  to  the  area  A  bounded  by  a  pair  of 
isochrons  and  a  pair  of  beam  boundaries.  We 


wish,  instead,  to  produce  a  scattering  strength  map 
in  terms  of  a  global  coordinate  system  and  thus 
we  will  use  Eq.  (1)  to  reexpress  the  reverberation 
intensity  in  terms  of  contributions  from  individual 
bathymetry  elements.  To  this  end,  we  again  use 
the  factored  form  of  the  scattering  function,  Eq. 
(4),  and  assume  that  all  quantities,  except  g*, 
associated  with  a  small  patch,  k,  of  the  bathymet¬ 
ric  grid  are  constant.  The  scattering  response  of 
patch  k  can  then  be  approximated  as 

dA 

«‘)  =  Y.  SW =  Y»  m,  G,(f)  (8) 

where,  for  notational  convenience,  we  have  incor¬ 
porated  the  beamformer  factor  b(<pb,(pa)  into  the  gen¬ 
eral  gain  factor  7k-  dA^/dt  is  the  rate  at  which  the 
patch  area  is  swept  out  by  the  isochrons;  it  is  zero 
for  t  <  Tjnin.k  and  t  >  Tmax.k,  the  minimum  and  maxi¬ 
mum  travel  times  associated  with  the  patch.  Convo¬ 
lution  of  g*  with  the  sweep  rate  function  acts  to 
smooth  g*,  suggesting,  once  again,  that  the  knowl¬ 
edge  of  the  details  of  the  scattering  is  not  critical 
provided  that  estimates  are  averages  over  large 
enough  areas,  a  necessary  condition  for  which  is 
that  Tjuax  V;  -  Tniin.k  »  W. 

The  total  scattered  intensity  of  a  single  beam  is 
the  sum  of  the  individual  patches  responses,  i.e. 

I(t)  =  Xlk(t)  =  EykmkGk(t)  (9) 

The  summation  implicitly  includes  the  left-right 
ambiguity  of  the  receiving  array  through  the  patch 
weighting  yk»  since  the  beamformer  factor  is  non¬ 
zero  only  along  the  conjugate  beam  directions.  As 
before  Eq.  (9)  could  be  explicitly  averaged  by  inte¬ 
gration  over  a  time  intervals  Ta  which  would  serve 
to  stabilize  the  individual  intensity  values  and  re¬ 
duce  the  size  of  the  subsequent  matrix  equations.  If 
the  scattering  function  is  taken  as  a  8-function,  av¬ 
eraging  yields 

pl{t)Ta= 'LykmkAkfk{t)  (10) 

where  A,,  is  the  area  of  the  kth  patch  and  fk  is  the 
fraction  of  the  area  of  the  kth  patch  ensonified  dur¬ 
ing  the  averaging  interval.  If  the  fractional  contri¬ 
butions  are  rounded  to  0  or  1 ,  then  Eq.  ( 10)  is  equiva¬ 
lent  to  Eq.  (5)  of  Ref.  (2).  In  the  latter,  a  synthetic 
test  case  was  considered  with  the  source  wavetrain 
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taken  as  2s  long  cw  pulse  and  the  averaging  time 
taken  equal  to  the  pulse  length,  resulting  in  mul¬ 
tiple  seafloor  patches  contributing  to  a  single  aver¬ 
aged  intensity  measurement.  In  practical  applica¬ 
tions,  though,  it  would  be  difficult  to  justify  the  use 
of  a  S-function  when  the  duration  of  the  source  pulse 
is  long  relative  to  the  size  of  the  individual  patches, 
since  implicitly,  as  noted  above,  it  ignores  the  con¬ 
volutional  nature  of  the  reverberation  process.  A 
scattering  function  whose  duration  matched  that  of 
Ae  cw  pulse  would  be  more  appropriate,  but  even 
so  the  structure  of  the  function  would  potentially  be 
important  for  good  scattering  estimates. 

A  direct  discretization  of  Eq.  (9)  for  time 
samples  att  =  jAt  is 

Ij  =  ?Ijk  =  ?GjkYkmk  (H) 

Individual  beam  responses  from  multiple  pings 
can  be  combined  to  form  a  single,  large  matrix  equa¬ 
tion  relating  reverberation  to  the  scattering  coeffi¬ 
cients  of  the  patches.  For  the  datasets  considered 
below,  the  number  of  scattering  coefficients  would 
be  -30, 000,  and  the  number  of  intensity  samples 
-20,000,  a  value  which  increases  proportionately  if 
multiple  pings  are  combined.  However,  the  matrix 
is  sparse  since  K(j),  the  number  of  patches  illumi¬ 
nated  at  time  sample,  j,  is  on  the  of  order  10  for  a 
single  ping. 

The  matrix  equation  could  be  solved  as  a  least 
squares  problem  with  the  addition  of  regularization 
constraints  on  the  scattering  coefficients,  using  any 
suitable  sparse  matrix  solver  such  as  conjugate  gra¬ 
dients10  or  LSQR1 1 ,  an  approach  we  investigate  else¬ 
where5.  Here  we  choose  to  solve  the  inverse  prob¬ 
lem  approximately  using  a  form  of  iterative 
backprojection.  The  advantage  of  backprojection  is 
that  it  is  simple,  fast  and  robust,  while  at  the  same 
time  capable  of  providing  a  fair  degree  of  ambigu¬ 
ity  resolution;  we  use  it  here  to  check  the  ping-to- 
ping  consistency  of  the  reverberation  data.  Although, 
backprojection  does  not  have  the  potential  resolu¬ 
tion  of  more  complete  inversion  methods,  especially 
for  multiple  pings,  these  methods  will  not  perform 
optimally  unless  the  data  are  essentially  consistent 
and  we  can  be  assured  that  there  are  no  unresolved 
problems  with  ship  position  and  array  orientation4. 

For  a  single  beam,  the  intensity  contribution  as¬ 
signed  to  a  patch  at  time  sample  j  by  the 


backprojection  is  proportional  to  the  weighting  fac¬ 
tor  for  that  patch 

Ijk  =  wkIj  (12) 

where 


Tk-  Gjk 


k^(i)Yk'Gik 


(13) 


The  backprojection  thus  incorporates  a  degree 
of  ambiguity  resolution,  assigning  intensity  based 
on  such  factors  as  transmission  loss  to  the  ensonified 
patches  and  also  on  the  basis  of  the  area  swept  out  - 
a  face  pointing  more  directly  into  the  acoustic  beam 
will  attract  proportionately  higher  assigned  inten¬ 
sity. 

The  result  of  the  backprojection  is  an  estimate 
of  the  scattering  intensity  function  Ik(t)  associated 
with  each  patch  k.  An  estimate  of  the  scattering  co¬ 
efficient,  mk,  is  found  by  integrating  the  scattering 
intensity,  Eq.  (9),  and  using  the  fact  that  the  integral 
of  g*(t)  is  unity  and  thus  the  integral  of  G(t)  is  Afc, 
the  area  of  patch  k.  The  discrete  time  version  of  the 
estimate  is 


pjk 

TkAk 


or  substituting  for  yk 


(14) 


1  [t„yXr'Zh 
Y W 


(15) 


This  latter  form  shows  explicitly  that  the  scattering 
coefficient  estimate  is,  as  it  should  be,  the  ratio  of 
the  scattered  intensity  to  the  incident  intensity,  scaled 
by  the  area. 

The  summation  over  the  individual  intensity 
contributions  helps  to  stabilize  the  scattering  coef¬ 
ficient  estimate  for  each  patch.  A  rough  estimate  of 
its  expected  variance  can  be  obtained  by  consider¬ 
ing  the  variance  of  the  contributing  intensity  val¬ 
ues,  assuming  the  wk  are  approximately  constant. 
The  time  bandwidth  product  for  the  LFM  pulse  and 
200  m  patches  is  approximately  15,  thus  the  expected 
variance  of  the  sum  is  approximately  1  dB’2,4. 

If  a  patch  is  illuminated  by  more  than  one  beam, 
then  the  scattering  coefficient  estimate  is  modified 


A.  J.  Harding,  M.  H.  Hedlin  and  J.A.  Orcutt:  Migration  of  backscatter  8 


appropriately  to  become  a  weighted  sum  of  the  in¬ 
dividual  beam  contributions. 


m* 


Ak 


(16) 


where  Ji^  is  the  fractional  area  of  patch  k  illumi¬ 
nated  by  beam  n  and  y^  is  the  appropriate  gain  fac¬ 
tor  for  the  beam. 

In  the  examples  examined  later,  we  have  em¬ 
ployed  an  iterative  form  of  the  backprojection  in 
which  the  weighting  function  in  Eq.  (12)  is  modi¬ 
fied  to 


wk  = 


Yk-  rak-Gkj 

k2o)Ykfik.Gkj 


(17) 


where  mk  is  the  scattering  coefficient  estimate  from 
the  previous  iteration.  The  numerator  is  now  the 
predicted  scattering  intensity  for  each  patch,  and  the 
effect  of  the  backprojection,  when  there  is  a  mis¬ 
match  between  the  predicted  and  recorded  intensi¬ 
ties,  is  to  modify  each  of  the  assigned  intensities  by 
an  equal  decibel  increment. 

The  backprojection  method  has  similarities  to 
the  "environmental  symmetry  breaking"  approach 
of  Ref.  MB,  particularly  on  the  first  iteration  when 
the  factor  controlling  the  assignment  of  energy  is 
usually  transmission  loss.  The  difference  between 
the  two  is  that  Ref.  MB  require  the  difference  in 
transmission  loss  to  reach  a  threshold  before  assign¬ 
ing  all  energy  to  one  side,  while  here  assignment  of 
energy  changes  continuously  as  a  function  of  the 
weighting  factors. 

The  above  formulation  is  essentially  a  local  one, 
relating  the  reverberation  signal,  sections  of  elevated 
intensity,  to  scattering  produced  by  direct  paths  to 
and  from  the  seafloor,  and  treats  non-local,  multiple 
scattering  as  part  of  the  reverberation  noise.  In  prin¬ 
ciple,  it  would  be  possible  to  incorporate  multiple 
scattering  via  suitable  paths.  However,  this  approach 
would  itself  require  an  a-priori  parameterization  of 
the  forward  scattering  and  would  quickly  become 
unwieldy  if  many  paths  were  included.  Similarly, 
long-range,  sub-seafloor  refraction  paths  are  as¬ 
sumed  not  to  contribute  significantly  to  the  rever¬ 
beration  signal,  although  local  volume  scattering  can 


'  be  considered  to  be  included  as  part  of  the  scatter¬ 
ing  function  g.  The  validity  of  the  local  scattering 
assumption  can  be  examined,  a-posteriori,  by  de¬ 
termining  whether  there  are  any  regions  of  high  scat¬ 
tering  coefficient  not  associated  with  identifiable 
bathymetric  features. 

III.  IMPLEMENTATION  DETAILS 

In  implementing  the  migration  algorithm  for  the 
SRP  data,  we  used  as  the  grid  for  the  scattering  co¬ 
efficient  map  the  200  m  by  200  m  swath  bathym¬ 
etry  grid  from  the  Hydrosweep  survey.  This  grid  size 
is  smaller  than  the  nominal  cross  track  resolution  of 
the  beamformed  data,  which  is  ~  600  m  at  broad¬ 
side  at  1/2  CZ,  but,  expressed  as  an  ensonification 
duration,  270-380  ms,  it  is  considerably  greater  than 
the  18  ms  duration  of  the  compressed  LFM  pulse. 
For  such  a  broadband  signal,  the  compressed  pulse 
width  yields  an  optimistic  estimate  of  range  resolu¬ 
tion,  since  the  actual  resolution  would  depend 
strongly  on  the  duration  of  the  scattering  response, 
including  the  surface  reflection  delay  at  the  source/ 
receiver.  Whether  the  200  m  grid  is  large  with  re¬ 
spect  to  the  scattering  function  duration  can  be 
checked  a-posteriori  by  assessing  the  stability  of  the 
scattering  strength  maps. 

The  quantities  needed  for  backprojection  of  a 
given  ping  are  calculated  in  two  stages,  first  relevant 
propagation  quantities  such  as  travel  time  and  ray 
angles  are  found  at  the  vertices  of  the  grid,  then  the 
sweep  function  G(t),  is  estimated  for  each  200  by 
200  m  grid  element,  which  we  refer  to  here  as  a 
patch.  The  spatial  homogeneity  of  the  sound  speed 
profile  permits  a  considerable  reduction  computa¬ 
tional  effort:  propagation  quantities  are  not  found 
directly  by  three-dimensional  ray  tracing  but  are  in¬ 
stead  interpolated  from  fields  calculated  on  regular 
two-dimensional,  range-depth,  grids.  This  simplifi¬ 
cation  also  makes  it  computationally  feasible  to  use 
PE  calculations,  rather  than  less  accurate  ray  theory, 
to  estimate  broadband  transmission  losses.  Two  sets 
of  PE  results  were  calculated,  one  with  a  ten  ele¬ 
ment  source  array  for  the  forward  direction,  the  other 
with  a  single  element  at  the  HLA  depth  for  the  re¬ 
verse  direction,  assuming  reciprocity.  The  single  fre¬ 
quency  results  were  then  integrated  across  the  source 
band  to  obtain  the  broadband  estimates. 

The  properties  of  the  individual  patches  are  cal- 
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culated  from  the  vertex  values.  The  patches  are  as¬ 
sumed  planar,  and  their  orientation  is  found  by  a 
least  squares  fit  to  the  vertex  depths.  Similarly  the 
isochron  surfaces  of  two-way  travel  time  are  as¬ 
sumed  to  be  locally  planar  and  are  found  by  a  least 
squares  fit  to  the  travel  times  at  the  vertices.  In  the 
following  examples  we  take  scattering  function,  g(t) 
to  be  a  5  function,  and  thus  the  sweep  function  G(t) 
for  a  patch  is  comprised  of  a  series  of  straight  line 
segments  with  vertices  corresponding  to  times  when 
the  wavefront  passes  a  patch  vertex  (Fig.  5).  From 
Fig.  5,  it  is  evident  that  using  a  more  complex  short 
duration  function  rather  than  a  8-function,  for  ex¬ 
ample  a  18  ms  boxcar,  would  have  only  a  minor 
effect  on  the  scattering  coefficient  estimates  after 
averaging. 

Geometric  shadowing  calculations  were  per¬ 
formed  prior  to  migration  using  a  mean  bathymet¬ 
ric  profile  along  each  beam  direction.  A  small  tran¬ 
sition  region  of  300  m  was  added  to  the  edge  of 
each  ensonified  region  in  order  to  help  prevent  small 
numerical  artifacts  or  local  roughness  from  produc¬ 
ing  shadow  zones.  Thus  individual  patches  that  point 
away  from  the  local  ray  direction  were  considered 
ensonified,  and  it  was  left  to  the  migration  to  deter¬ 
mine  whether  they  produced  small  backscatter.  Dur¬ 
ing  backprojection,  the  scattering  strengths  of 
patches  in  the  shadow  zones  were  adjusted  at  each 
iteration  in  order  to  produce  an  expected 
backscattered  intensity  equal  to  the  10th  percentile 
of  the  ensonified  patches.  The  purpose  of  this  was 
to  ensure  that  the  preponderance  of  energy  would 
be  assigned  to  an  ensonified  area  if  the  conjugate 
side  was  in  shadow,  while  at  the  same  time  permit¬ 
ting  energy  assignment  if  both  sides  were  in  shadow. 
In  this  way  it  was  possible  to  keep  track  of  unex¬ 
pected  scattering  sites  in  the  reverberation  map. 

IV.  MIGRATION  EXAMPLE  FROM  B’ 

We  take  as  an  initial  migration  example,  ping 
1313  of  segment  s415,  for  which  the  RN  Cory 
Chouest  was  located  just  north  and  west  of  the  cen¬ 
ter  of  the  track  star  near  ridge  B\  Fig.  1,  and  had  a 
HLA  heading  of  169°,  Table  I.  At  a  1/2  CZ  on  the 
starboard  side  is  B’  with  its  backfacing  scarps  ori¬ 
ented  almost  perpendicular  to  the  beam  directions, 
while  at  the  same  range  on  the  port  side,  a  series  of 
N25°E  trending  ridges  point  almost  directly  into  the 


FIG.  5.  Example  sweep  area  function  for  a  seafloor  patch. 
The  sweep  function  consists  of  linear  segments  with  joints 
corresponding  to  the  passage  of  the  beam  wavefront  (dash 
lines)  past  the  patch  vertices.  For  the  calculation  both  the  patch 
and  wavefront  are  assumed  locally  planar.  (Patch  area  is  nor¬ 
malized  to  unit  area). 

beam  directions  (Fig.  1).  The  reverberation  returns 
from  these  ridges  and  the  scarp  of  B’  interfere  to 
produce  a  set  of  crossing,  high  amplitude  events  in 
the  data  between  30  and  45  s,  Fig.  6.  When  we  ex¬ 
amine  the  two-way  transmission  loss  to  the  conju¬ 
gate  areas,  we  find  that  only  the  port  side  ridges  and 
the  upper  part  of  B’  are  shallow  enough  to  project 
into  the  center  of  the  acoustics  beam,  although  a 
broad  platform  on  the  port  side  is  also  well 
ensonified,  Fig.  7. 

Fig.  8  displays  the  reverberation  map  obtained 
after  the  first  iteration  of  the  migration  using  data 
from  20s  out  to  times  corresponding  to  the  1/2  CZ 
turning  ranges.  The  reverberation  map  is  an  inter¬ 
mediate  backprojection  output  and  is  calculated  from 
the  sum  of  the  backscattered  energy  associated  with 
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Pinq  1313:  Time  series 
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FIG.  6.  (a)  Reverberation  data  for  ping  1333,  broadband  LFM  transmission  of  s415.  The  displayed  data  spans  the  time  interval 
of  reverberation  returns  from  1/2  CZ  and  beam  numbers  15-120,  which  avoids  endfire  (Beam  nos.  64/65  are  the  broadside 
beams).  To  facilitate  comparision  with  the  reverberation  maps,  the  data  has  been  averaged  over  an  interval  of  270  ms  and 
corrected  for  the  average  number  of  ensonified  patches,  a  correction  of  -1 1 .1  dB.  The  prominent  returns  between  30  &  45s  are 
a  combination  of  returns  from  scarps  on  the  front  face  of  B’  on  the  starboard  side  of  the  array  and  ridges  on  the  port  side  that 
point  almost  directly  into  the  beams,  (b)  &  (c)  estimated  port  and  starboard  lime  series  after  iterative  migration.  The  split  is 
plausible  and  for  the  most  part  clean,  although  some  residual  energy  appears  to  be  misplaced. 
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FIG.  7.  Two  way  transmission  loss  for  ping  1313  projected  onto  the  port  and  starboard  side  bathymetry.  Strongly  ensonified 
features  include  the  top  of  B'  on  the  starboard  side,  and  the  high  ridges  pointing  approximately  into  the  beams  on  the  port  side. 
Current  figure  does  not  take  into  account  geometric  shadowing. 
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FIG.  8.  Reverberation  map  from  the  initial  migration  of  ping  1313.  The  left-right  ambiguity  resolution  is  based  primarily  on 
ensonification  level,  although  some  dependence  on  local  incidence  angle  is  included  via  the  sweep  rale  function.  High  amplitude 
returns  due  to  local  slope  such  as  those  from  the  scarp  slopes  of  B’  are  not  fully  resolved  and  also  appear  in  the  conjugate  image. 


A.  J.  Harding,  M.  H.  Hedlin  and  J.A.  Orcutt:  Migration  of  hackscatter  12 


each  patch  of  the  grid,  Eq.  (12).  The  reverberation 
map  is  a  stable  estimate  in  the  sense  that 
backprojection  is  guaranteed  to  be  intensity  preserv¬ 
ing,  and  the  partition  of  energy  between  simulta¬ 
neously  ensonified  patches  is  affected  only  by  rela¬ 
tive  not  absolute  errors  in,  for  example,  transmis- 
v':*  sion  loss.  From  Fig.  8  we  can  see  that  the  high  am¬ 
plitude  backscatter  events  are  predominantly  located 
within  areas  that  were  strongly  ensonified.  However, 
a  notable  exception  is  the  strong  scattering  associ¬ 
ated  with  a  line  of  small  scarps  on  B'  in  the  region 
(2955-2970,  205-210)  UTM  km,  which  lies  well 
outside  the  main  acoustic  beam.  Furthermore,  it  is 
evident  that  within  the  strongly  ensonified  regions 
the  reverberation  events  tend  to  be  localized  along 
bathymetric  features  with  high  slopes,  such  as  the 
scarps  covering  the  upper  parts  of  B\  This  impres¬ 
sion  is  supported  by  Fig.  9,  which  shows  that  above 
a  cutoff  of  about  10°  the  mean  scattering  strength 
increases  with  local  grazing  angle,  a  dependence  that 
arises  even  though  the  a-priori  assumption  is  that 
the  scattering  coefficient  is  constant.  Local  slope  in¬ 
fluences  the  initi  al  partition  coefficients  only  through 
its  effect  on  the  sweep  rate  function  for  each  patch. 

Although,  the  initial  migration  produces  a  fair 
'  degree  of  left-right  ambiguity  resolution,  it  is  clearly 
not  perfect  if  local  slope  is  the  sole  or  primary  de- 
terminant  of  backscatter  strength.  For  example,  en¬ 
ergy  that  on  the  starboard  side  is  associated  with 
the  scarps  of  B’,  also  appears  on  the  port  side  where 
it  cuts  across  the  mostly  shallow  slope  bathymetric 
fabric.  We  can  test  the  degree  to  which  slope  ac¬ 
counts  for  strong  backscatter  by  incorporating  the 
estimated  angular  dependence  of  the  scattering  co¬ 
efficient  into  the  partition  coefficients  for  the 
backprojection,  Eq.  (17).  All  events  associated  with 
high  slope  features  will  migrate  solely  to  the  appro¬ 
priate  left  or  right  image,  only  anomalous  events 
will  remain  ambiguous  and  appear  on  both  images. 
For  this  dataset,  the  migration  effectively  converged 
after  three  iterations,  with  only  a  barely  perceptible 
change  in  the  mean  backscatter  coefficient  between 
iterations  3  and  4,  Fig.  9.  The  perceived  ambiguity 
in  the  backscattered  energy  images  is  significantly 
reduced  and  only  a  few  small  isolated  events  ap¬ 
pear  split  between  the  two  images,  Fig.  10.  The  ef¬ 
fectiveness  of  the  ambiguity  resolution  is  more 
readily  assessed  in  the  original  time-beam  coordi- 


FIG.  9.  Estimate  of  backscatter  strength  as  a  function  of  graz¬ 
ing  angle  for  the  4  migration  iterations  of  Ping  1313.  Iteration 
1-  dotted,  2-  dash-dot,  3-dashed,  4-solid.  For  iteration  1  there 
is  no-apriori  assumption  of  increased  backscatter  with  graz¬ 
ing  angle,  yet  the  mean  backscatter  strength  shows  a  correla¬ 
tion  with  grazing  angle  down  to  about  10°.  For  subsequent 
iterations  the  mean  backscatter  of  the  previous  iteration  is  used 
and  convergence  occurs  in  3-4  iterations. 

nates,  Fig.  6.  For  the  most  part  this  view  confirms 
the  impression  that  the  energy  is  cleanly  split  be¬ 
tween  the  left  and  right  side.  However,  energy  near 
the  edges  of  some  prominent  reverberation  events 
leaks  through  to  the  other  side,  suggesting  that  there 
may  be  a  slight  error  in  the  array  orientation  or  that 
the  migrated  image  could  be  improved  by  incorpo¬ 
rating  some  cross  talk  between  beams  in  the 
beamformer  factor  b,  Eq.  (4). 

The  scattering  strength  map  for  the  4th  iteration 
is  displayed  in  Fig.  1 1.  As  expected  the  elongated 
scarp  slopes  on  B’  are  highlighted  in  the  image,  but 
compared  to  the  energy  image,  the  prominence  of 
the  port  ridges  is  reduced,  their  previous  prominence 
being  due  solely  to  strong  ensonification.  As  was  to 
be  expected  from  the  mean  behavior,  there  is  an 
overall  correspondence  between  the  grazing  angle 
and  scattering  strength  with  the  normalization  em¬ 
phasizing  the  scattering  strength  of  the  lower  slopes. 
Also,  predictably,  the  scattering  strength  map  is 
noisier  than  the  reverberation  one,  since  now  global 
as  well  as  local  errors  in  estimated  quantities  are 
important.  For  example,  any  mismatch  between  the 
predicted  and  actual  pattern  of  transmission  loss  will 
be  reflected  in  the  scattering  strength  map. 
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FIG.  10.  Reverberation  map  of  ping  1313  after  four  iterations  .  Compared  to  iteration  1,  the  left-right  ambiguity  resolution  is 
significantly  improved,  with  a  much  weaker  ghosting  of  the  B'  scarps  in  the  port  side  image. 
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FIG.  ll.Scattering  strength  maps  derived  from  the  reverbearation  map,  Fig  10.,  by  correcting  for  transmission  loss,  area,  etc. 
The  previous  prominence  of  the  top  of  B'  and  the  port  side  ridges  is  revealed  as  primarily  a  consequence  of  high  ensonification 
levels,  while  the  B*  scarp  slopes  are  identified  clearly  as  regions  of  high  backscauer. 
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V.  COMPOSITE  IMAGES  OF  B’ 

We  can  assess  the  validity  of  the  single  ping 
ambiguity  resolution  and  simultaneously  gain  in¬ 
sight  into  the  general  quality  of  the  data  by  compar¬ 
ing  the  migration  results  from  multiple 
ensonifications  of  the  same  target  with  nearly  the 
same  ship  location  but  different  headings.  The  stan¬ 
dard  means  of  obtaining  multiple  looks  at  a  target  is 
to  use  reverberation  records  from  crossing  paths  and 
this  we  do  by  comparing  images  from  runs  5a  and  6 
of  the  MAE  which  cross  at  an  angle  of  approxi¬ 
mately  25°.  But  we  can  also  gain  useful  informa¬ 
tion  from  multiple  pings  from  a  single  run,  taking 
advantage  of  small  changes  in  the  heading  of  the 
receiver  array  between  pings.  If  a  reverberation  event 
is  migrated  to  its  true  origin  on  the  seafloor,  then  its 
location  will  not  change  when  the  array  heading 
changes  by  0.  Conversely  a  false  scatterer  that  re¬ 
sulted  from  migrating  energy  to  the  wrong  side  of 
the  array  will  rotate  through  20  since  it  is  located  at 
the  mirror  image  of  true  scatterer  in  the  array,  Fig. 
12.  Thus  the  incorrectly  migrated  energy  will  move 
in  a  predictable  fashion  in  step  with  the  changes  in 
array  orientation.  This  movement  can  most  readily 
be  appreciated  by  combining  the  migration  results 
from  successive  pings  into  an  animation  sequence. 
For  the  5  consecutive  LFM  pings  from  run  5a  con¬ 
sidered  here,  the  maximum  difference  in  array  head¬ 
ing  is  3°  which  is  sufficient  to  move  false  scatterers 
across  3-4  beams  in  the  migrated  image. 

We  have  investigated  the  consistency  of  the  mi¬ 
grated  images  using  5  consecutive  broadband  LFM 
pings  from  both  run  5a  and  run  6  spanning  their 
intersection  at  26°  36’  N,  -47°  47'  W  ( (2946, 222), 
Fig.  1.).  Fig.  13  displays  the  scattering  strength  maps 
for  the  northern  foreslopes  of  B’  for  9  of  the  pings. 
Perhaps  the  most  striking  aspect  of  the  images  is 
that  relatively  narrow  scarp  slopes  appear  consis¬ 
tently,  with  only  relatively  minor  ping-to-ping  varia¬ 
tion  in  strength  and  location.  The  images  from  run 
6  (pings  1683-1696)  are  cleaner  than  the  ones  from 
run  5a  (pings  1313-1336)  due  to  the  fact  that  the 
conjugate  area  on  run  6  is  the  source  of  much  less 
backscatter.  For  run  5a,  the  lower  portion  of  the 
images  below  2950  includes  residual  mislocated 
energy  from  the  port  side  ridges.  A  port  side  origin 
is  supported  by  the  fact  that  it  does  not  appear  in  the 
images  from  run  6  and  that  its  location  moves  be¬ 
tween  the  successive  run  5a  images.  Similarly  en- 


FIG.  12.  With  incomplete  ambiguity  resolution,  energy  from 
a  scatterer  will  also  be  migrated  to  a  false  scatterer  location; 
that  is,  the  mirror  image  in  the  line  array  of  the  true  location. 
These  locations  will  remain  paired  for  different  pings  (black 
circles)  if  the  ship  heading  remains  constant,  but  the  location 
of  the  false  scatterer  will  rotate  through  20  when  the  array 
heading  changes  by  0. 

ergy  that  appears  in  the  run  5a  images  near  the  cen¬ 
ter  of  a  small  basin  at  (2960,202)  can  also  be  attrib¬ 
uted  to  a  partial  failing  of  the  ambiguity  resolution. 

Averaging  the  individual  images  effectively  sup¬ 
presses  false  scatterers  attributable  to  energy  leak¬ 
age  and  emphasizes  the  consistency  of  high  back¬ 
scatter  from  the  narrow  ridges,  Fig.  14.  All  the  high 
backscatter  regions  correspond  to  scarp  features  that 
are  ensonified  at  high  local  grazing  angles.  How¬ 
ever,  the  converse  is  not  true,  there  are  features,  most 
notably  the  front  edge  of  a  small  platform  centered 
(2955, 198)  that  although  illuminated  at  high  angle 
does  not  produce  high  backscatter.  Spots  of  high 
backscatter  do  appear  along  the  front  edge  of  the 
platform  in  the  individual  migration  images,  per¬ 
haps  indicative  of  backscatter  speckle  that  fluctu¬ 
ates  in  response  to  changes  in  the  ensonification. 
However,  since  the  strongest  backscatter  events  ap¬ 
pear  in  the  noisier  run  5a  images  it  is  more  likely 
that  the  high  ensonification  angles  succeed  in  cap¬ 
turing  port  side  energy  during  migration. 

For  each  of  the  10  pings,  we  have  used  the  mi¬ 
gration  results  to  estimate  the  variation  of  mean 
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Fig.  13.  Scattering  coefficient  maps  of  a  subarea  of  the  B'  abyssal  hill  for  5  consecutive  broadband  LFM  pings  of  run  5a,  1313 
through  1336,  and  4  consecutive  pings  for  run  6, 1683-1696.  Most  of  the  scarp  slopes  in  the  area  show  up  consistently  as  regions 
of  high  backscatler.  A  notable  exception  is  the  edge  of  a  small  platform  near  (2955,200)  which  is  qnly  intermittently  highlighted 
possibly  only  as  the  result  of  being  mistakenly  assigned  energy  from  the  conjugate  area. 
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Fig.  14.  Comparison  of  the  mean  backscatter  strength,  left  panel,  found  by  stacking  the  10  migrated  pings,  with  the  mean 
illumination  angle,  right  panel,  expressed  in  terms  of  the  cosine  of  the  local  grazing  angle.  For  the  most  part,  there  is  a  close 
correspodance  between  high  backscatter  targets  and  local  grazing  angle:  all  prominent  backscatter  targets  correspond  to  scarp 
features  ensonified  at  high  grazing  angles.  However,  there  are  some  features,  most  notably  the  front  slope  of  a  small  platform 
centered  (2955, 198)  that  although  illuminated  at  high  angle  does  not  produce  a  high  backscatter. 


backscatter  strength  with  grazing  angle.  Fig.  1 5.  The 
backscatter  curve  was  found  using  the  Loess  algo¬ 
rithm13,  assuming  that  the  local  grazing  angles  were 
correct.  The  Loess  algorithm  finds  a  point  on  the 
backscatter  curve  for  a  given  grazing  angle  by  fit¬ 
ting  a  line  to  a  subset  of  the  data  centered  on  the 
grazing  angle  using  weighted  least  squares.  Formal 
confidence  intervals  for  the  mean  were  estimated 
from  200  bootstrap  samples14  for  each  ping,  and  were 
typically  1-2  dB  at  the  95  %  confidence  level,  in¬ 
creasing  at  high  grazing  angles.  Over  the  interval 
from  10°  to  45°,  the  mean  scattering  strength  curves 
for  run  5a  are  roughly  comparable  to  a  Lambert’s 
law  type  behavior  with  a  Mackenzie  parameter  of  - 
27  dB,  although  the  curves  are  distinct  based  on  the 
formal  confidence  intervals.  The  backscatter  esti¬ 
mates  for  a  given  grazing  angle  are  distributed  ap¬ 


proximately  normally  with  a  standard  deviation  of 
4-5  dB,  Fig.  16. 

The  split  is  probably  a  consequence  of  differ¬ 
ences  in  the  distribution  of  the  relatively  small  num¬ 
ber  of  patches  that  are  ensonified  at  high  grazing 
angles.  For  run  5a,  73%  of  the  700-900  patches 
ensonified  at  angles  greater  than  25°  (approximately 
3%  of  the  total)  come  from  the  B’  scarps,  while  for 
run  6,  this  percentage  drops  to  63%  with  a  greater 
number  of  high  angle  targets  on  the  starboard  side 
that  are  small  and  relatively  isolated.  When  the  back¬ 
scatter  curves  are  split  into  port  and  starboard  con¬ 
tributions,  the  gap  between  the  starboard  side  curves 
that  includes  B’  is  reduced,  although  not  eliminated. 
Moreover,  the  port  side  targets  display  a  weaker 
dependence  on  angle.  Fig.  17,  and  are  consistent 
over  the  two  runs.  The  difference  in  scattering 
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FIG.  15.  Estimate  of  scattering  strength  vs.  grazing  angle  for 
pings  from  run  5a  (solid)  and  run  6  (dashed)  with  bootstrap 
estimates  of  ±2s  errors  in  the  mean.  For  reference,  Lamberts 
law  curve  using  a  Mackenzie  parameter  of  -27  dB  is  included, 
dotted  curve. 

strength  may  indeed  reflect  an  intrinsic  difference 
in  the  backscatter  level  between  the  areas.  However, 
it  may  also  be  a  consequence  of  processing  artifacts 
and  of  the  limitations  of  the  bathymetry  database. 
The  size  of  the  port  side  targets  tend  to  be  near  or 
below  the  angular  resolution  of  the  system,  and  thus 
the  intrinsic  averaging  of  the  beamforming  coupled 
with  any  small  errors  in  the  heading  calculations 
will  tend  to  mute  the  influence  of  high  angle  patches. 
Also  the  number  of  high  angle  targets  is  small 
enough  that  imperfections  in  the  bathymetry  data¬ 
base  such  as  small  noise  spikes  that  produce  erro¬ 
neous  high  angle  targets  would  bias  the  scattering 
strength  estimates  downward. 

VI.  CONCLUSIONS 

We  have  developed  an  algorithm  for  migrating 
reverberation  data  that  is  based  upon  time-domain 
expressions  for  the  expected  intensity  of  the  rever¬ 
beration  field.  The  time  dependent  aspects  of  the 
field  are  parameterized  by  a  function  referred  to  as 
the  scattering  function,  which  includes  both  stochas¬ 
tic  components  dependent  on  the  details  of  the  sea¬ 
floor  scattering  and  deterministic  ones  due  to  the 
characteristics  of  the  acoustic  imaging  system.  We 
argue  that  the  details  of  the  scattering  function  are 
not  critical  provided  that  scattering  coefficients  are 
estimated  over  patches  of  the  seafloor  with  response 


FIG.  16.  Representative  histogram  of  scattering  strength 
relavite  to  mean  trend  from  ping  1 3 1 3  for  a  bin  extending  from 
21°-29°.  There  are  1 385  patches  in  the  bin.  The  distribution  is 
approximately  gaussian  with  a  =  4.4  dB 

times  that  are  long  compared  with  the  characteristic 
time  of  the  scattering  function,  or  equivalently  whose 
area  is  large  compared  to  the  correlation  scales  of 
the  scattering  process9.  This  approach  is  comparable 
to  reducing  speckle  in  optics  by  averaging  over  a 
scanning  aperture15.  Time  averaging  has  the  addi¬ 
tional  advantage  of  producing  stable  estimates  from 
individual  reverberation  records,  the  trade-off  is  that 
spatial  resolution  is  reduced.  However,  even  with 
averaging,  the  broadband  LFM  pulse  of  the  SRP 
experiments  is  capable  of  yielding  consistent  scat¬ 
tering  strength  estimates  at  the  resolution  of  the  main 
bathymetry  database. 

A  key  feature  of  the  algorithm  is  that  the  results 
are  computed  in  the  global  coordinate  system  of  the 
bathymetry  database  rather  than  the  beam-travel 
time  coordinates  of  individual  pings,  facilitating  the 
comparison  of  results  from  multiple  pings.  The  ba¬ 
sis  for  the  coordinate  conversion  is  the  precise, 
sample-by-sample,  tracking  of  the  acoustic 
wavefronts  across  the  individual  patches  and  the 
assumption  that  the  scattering  coefficient  is  constant 
within  individual  elements,  patches,  of  the  bathym¬ 
etry  database.  These  two  assumptions  allow  us  to 
reexpress  the  reverberation  response  as  a  large  sparse 
matrix  equation  with  the  scattering  coefficients  as 
the  unknowns.  The  matrix  formulation  enables  mul¬ 
tiple  records  to  be  combined  simply  into  a  single 
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FIG.  17.  Port  and  starboard  estimates  of  mean  scattering 
strength  for  the  10  pings.  Estimates  from  run  5a  are  solid  lines 
and  those  from  run  6  are  dashed.  The  estimates  from  the  star¬ 
board  side  that  ensonify  B’  show  a  consistently  stronger  de¬ 
pendence  on  grazing  angle  and  there  is  less  of  a  gap  between 
run  5a  &  6  estimates. 

linear  inversion  problem. 

In  this  paper  we  have  taken  a  processing  ori¬ 
ented  approach  taking  advantage  of  the  sparsity  to 
solve  the  matrix  equation  approximately  using  it¬ 
erative  backprojection  rather  than  attempting  to 
solve  the  equation  directly  as  part  of  a  linear  inver¬ 
sion.  We  describe  this  process  as  migration  since 
the  method  is  analogous  to  the  Kirchoff  migration 
procedure  of  seismic  processing,  the  principal  dif¬ 
ference  being  that  here  it  is  applied  to  incoherent 
scattering  rather  than  coherent  reflections.  It  has  the 
processing  virtue  of  producing  reverberation  and 
scattering  coefficient  maps  rapidly  thus  facilitating 
efficient  analysis  of  multiple  records.  The  reverbera¬ 
tion  maps,  which  are  produced  by  backprojection 
as  an  intermediate  output,  represent  an  exact  parti¬ 
tion  of  the  record  intensity  and  are  thus  guaranteed 


to  be  stable  and  allow  us  to  assess  directly  the  ef¬ 
fectiveness  of  the  left-right  ambiguity  resolution. 
The  performance  of  the  migration  algorithm  is 
illustrated  using  two  sets  of  broadband 
monostatic  reverberation  records  from  a  pair  of 
crossing  runs  of  the  SRP  experiment  on  the 
flanks  of  the  Mid- Atlantic  ridge.  At  a  1/2  CZ, 
iterative  migration  produces  well  resolved  maps 
of  reverberation  and  scattering  strength.  In  the 
reverberation  maps,  there  is  a  clear  correspon¬ 
dence  between  the  location  of  high  amplitude 
reverberation  events  and  bathymetric  features 
with  either  low  transmission  loss,  high  local 
grazing  angle  or  both.  In  particular,  narrow  scarps 
less  than  1  km  wide  are  highlighted  consistently. 
The  migration  also  produces  a  plausible  resolu¬ 
tion  of  the  left-right  ambiguity,  dividing  the 
majority  of  events  in  a  reverberation  record 
cleanly  between  the  two  ensonified  areas.  The 
migration  algorithm  can  be  regarded  as  an  envi¬ 
ronmental  symmetry  breaking  method  that  uses 
seafloor  induced  asymmetries  to  resolve  the 
inherent  left-right  ambiguity  of  the  receiving 
array.  Although  the  migration  images  are  by  no 
means  unique,  the  improvement  of  the  event 
separation  with  iteration  and  the  clear  correspon¬ 
dence  between  events  and  bathymetric  features 
provides  strong  a-posteriori  justification  for  the 
inclusion  of  mean  backscatter  strength  as  a  factor 
in  the  backprojection  weights. 

The  results  of  multiple  migrations  from  a  single 
run  and  from  crossing  runs  corroborates  the  over¬ 
all  effectiveness  of  the  left-right  ambiguity  resolu¬ 
tion  for  individual  records,  and  supports  the  idea 
that  scattering  strength  is  primarily  a  function  of 
grazing  angle.  The  majority  of  high  scattering 
strengths  areas,  primarily  associated  with  scarp 
slopes,  appear  consistently  in  all  migration  images. 
A  few  potentially  anomalous  areas  move  in  step 
with  changes  in  array  heading  within  a  single  run 
and  tend  to  be  absent  in  images  from  the  crossing 
run.  These  areas  are  identified  as  energy 
misassigned  to  the  wrong  side  of  the  array,  and  are 
attenuated  when  multiple  maps  are  stacked  together. 
Employing  multiple  heading  directions  from  cross¬ 
ing  runs  is  the  standard  means  of  reducing  left-right 
ambiguity.  The  present  results  demonstrate  that 
even  small  heading  changes  (e.g  2-3°)  from  a  single 
run  can  be  profitably  employed  to  reduce  ambigu- 
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ity.  The  consistency  of  the  migration  images  dem¬ 
onstrates  the  stability  of  the  acoustic  imaging  sys- 
.  tern  used  for  the  SRP  experiments  and  is  a  good 
indication  that  linear  inversion  could  be  applied  suc¬ 
cessfully  to  the  data. 

The  migration  results  indicate  that  transmission 
v-:;  loss  and  grazing  angle  are  the  dominant  factors  con¬ 
trolling  reverberation  response  in  the  study  area. 
However,  weaker  dependencies  of  scattering 
strength  on  spatial  and  azimuthal  variations  in  scat¬ 
tering  strength  cannot  be  ruled  out  by  the  present 
analysis.  A  dependence,  in  the  mean,  of  backscatter 
strength  on  grazing  angle  appears  after  the  first  it¬ 
eration  of  the  migration  when  the  backprojection 
weights  contain  no  bias  towards  patches  with  high 
grazing  angles.  We  thus  conclude,  as  have  previous 
investigators,  that  the  correlation,  in  the  mean,  be- 
-  tween  backscatter  strength  and  grazing  angle  is  a 
robust  feature  of  the  data. 

Above  a  noise  floor  of  about  5°,  the  estimated 
dependence  of  the  mean  scattering  strength  on  graz¬ 
ing  angles  is  approximately  equal  to,  but  formally 
.  statistically  distinct  from,  that  of  Lambert’s  law  with 
a  Mackenzie  parameter  of  -27  dB.  Although  the 
general  magnitude  of  the  increase  in  backscatter 
: strength  appears  robust,  the  detailed  dependence  on 
grazing  angle  must  be  treated  with  some  caution  as 
the  accuracy  of  the  slopes  derived  from  the  swath 
bathymetry  is  itself  questionable,  especially  when 
the  high  slopes  are  associated  with  smaller  features 
near  the  resolution  limit  of  the  system16. 
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